 function f     = jacobLogit(theta,delta,Data)

% jacobLogit - This function is calculate the derivative of the mean value
% delta wrt the parameter using the implicit function theorem (see Nevo
% http://faculty.wcas.northwestern.edu/~ane686/supplements/Ras_guide_appendix.pdf
% 
% 
% jacobLogit(theta,delta,Data)
%
% Inputs:
%    input1 - Non linear parameters
%    input2 - Data
%
% Outputs:
%    Derivative of Delta wrt theta 
%
% Subfunctions: LogitShareCalculation; multinv; multiprod
%

% Author: Laura Grigolon
% August 2012;
 
% Unpack
prods           = Data.prods;
multixvuLRC     = Data.multixvuLRC;
multixvuLSEG    = Data.multixvuLSEG;
qweightrprods   = Data.qweightrprods;

% Market Shares
if numel(theta) == 1
    [~, sij]        = LogitShareCalculation(theta,delta,Data);
elseif numel(theta) == 2
    [~, sij]        = LogitShareCalculation2RC(theta,delta,Data);
end

% Jacobian
% remember that sij is prod x 1 x mkt x nodes

% derivative of shares wrt to delta
part1           = bsxfun(@times,sij,eye(prods));          % diagonal in multiple dimensions
sijtransp       = permute(sij,[2 1 3 4]);
part2           = multiprod(sij,sijtransp);
derShareDeltij  = part1 - part2;
derShareDelta   = sum((bsxfun(@times,derShareDeltij,qweightrprods)),4);

% derivative of shares with respect to random coefficient 1
if size((theta),1) == 1
    part1           = bsxfun(@times,multixvuLRC,sij);
    part2           = bsxfun(@times,sum(part1),sij);
    derShareRC1     = part1-part2;
    derShareTheta   = sum((bsxfun(@times,derShareRC1,qweightrprods)),4);
    derDeltaTheta   = multiprod((-multinv(derShareDelta)),derShareTheta);
    derDeltaTheta   = reshape(derDeltaTheta,[Data.nobs 1]); % 2D market shares
elseif size((theta),1) == 2
    % derivative of shares with respect to random coefficient on SEG
    part1RC1           = bsxfun(@times,multixvuLSEG,sij);
    part2RC1           = bsxfun(@times,sum(part1RC1),sij);
    derShareRC1        = part1RC1-part2RC1;
    derShareTheta1     = sum((bsxfun(@times,derShareRC1,qweightrprods)),4);
    % derivative of shares with respect to random coefficient on continuous
    % characteristic
    part1RC2           = bsxfun(@times,multixvuLRC,sij);
    part2RC2           = bsxfun(@times,sum(part1RC2),sij);
    derShareRC2        = part1RC2-part2RC2;
    derShareTheta2     = sum((bsxfun(@times,derShareRC2,qweightrprods)),4);
    
    derShareTheta      = [derShareTheta1 derShareTheta2];
    
    derDeltaTheta      = multiprod((-multinv(derShareDelta)),derShareTheta);
    derDeltaTheta      = [(reshape(derDeltaTheta(:,1,:),[Data.nobs 1])),(reshape(derDeltaTheta(:,2,:),[Data.nobs 1]))]; % 2D market shares
end

f = derDeltaTheta;

end